package modeling;


import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.File;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.text.DateFormat;
import java.text.DecimalFormat;
import java.text.NumberFormat;
import java.text.SimpleDateFormat;
import java.util.Date;
import java.util.Random;

import utilities.PDButility;


/** ------------------------------------------------------------------------------------------------------------------
 * start of MCMC5C class
 * Version 1.0
 * January 14th, 2011
 * ------------------------------------------------------------------------------------------------------------------  */
public class MCMC5C {


	/*  global variables  */
	private Date date;
	private DateFormat dateFormat = new SimpleDateFormat("yyyy/MM/dd\tHH:mm:ss");
	private final DateFormat dateFormatOutputFolder = new SimpleDateFormat("yyyy-MM-dd-HH-mm-ss");
	private final String outputFolderDate = dateFormatOutputFolder.format(new Date());
	private final String outputFolder = "MCMC5C/Output/" + outputFolderDate;
	private final String inputPath = "MCMC5C/Input/";
	private final String outputPath = "MCMC5C/Output/" + outputFolderDate + "/";
	private NumberFormat formatter = new DecimalFormat("0");
	private static PDButility pdb = new PDButility();
	private String filenameIFs = "no_name.txt";
	private String filenameFragments = "no_name.txt";
	private int numberOfPoints = 0;
	private boolean determinedNumberOfPoints = false;
	private int[] pointGenomicLocations;
	private int numberOfIFValues = 0;
	private int numberOfVarianceValues = 0;
	private double[][] IFvalues;
	private double[][] IFvariances;
	private double averageIF = 0;
	private double defaultVarianceValue = 0;
	private double[][] structure;
	private double[][] bestStructure;
	private double probabilityOfBestStructure;
	private long indexOfBestStructure;
	private long runLength = 1000000;
	private double numberOfStructuresToOutput;
	private double stepSize = 0.5;
	private String acceptanceFunctionType = "mcmc";			//  mcmc is the default value, optionally can use "gradient" as the acceptance function (simulate 5C3D method)
	private boolean acceptanceFunctionMCMC = true;
	private boolean acceptanceFunctionGradient = false;
	private long numberOfChangesAccepted = 0;
	private String[] argsGlobal;
	/*  end of global variables  */
	
	
	/*  start of getNumberOfPoints method  */
	private void getNumberOfPoints(String fileToCountNumberOfPointsFrom){
		String line = null;
		numberOfPoints = 0;
		try{
			BufferedReader br = new BufferedReader(new FileReader(inputPath + fileToCountNumberOfPointsFrom));			
			while((line = br.readLine()) != null){
				numberOfPoints++;
			}
			br.close();
		}catch (IOException e) {
			e.printStackTrace();
			System.out.println("Error: cannot find file with filename \"" + fileToCountNumberOfPointsFrom + "\"\n" +
			"System.exit(0)");
			System.exit(0);
		}
		determinedNumberOfPoints = true;
	}/*  end of getNumberOfPoints method  */

	
	/*  start of getFragmentsFromFile method  */
	private void getFragmentsFromFile(){
		//  make sure that the number of points has already been determined
		if(!determinedNumberOfPoints){
			System.out.println("Error: the number of points must be determined before reading the fragments as input.\nSystem.exit(0)");
			System.exit(0);
		}
		
		pointGenomicLocations = new int[numberOfPoints];

		String line = null;
		int lineCounter = 0;
		try{
			BufferedReader br = new BufferedReader(new FileReader(inputPath + filenameFragments));			
			while((line = br.readLine()) != null){
				String values[] = line.split("\t");
				if(values.length == 1){
					pointGenomicLocations[lineCounter] = Integer.parseInt(values[0]);
				}
				else{
					System.out.println("Error: there is not one entry in each row in the fragments file\n"
							+ "System.exit(0)");
					System.exit(0);
				}
				lineCounter++;
			}
			br.close();
		}catch (IOException e) {
			e.printStackTrace();
			System.out.println("Error: cannot find feature file with filename \"" + filenameFragments + "\"\n" +
			"System.exit(0)");
			System.exit(0);
		}
		if(lineCounter != numberOfPoints){
			System.out.println("Error: there are not the same number of fragments (" + lineCounter + ") as there are points (" + numberOfPoints + ")\n"
					+ "System.exit(0)");
			System.exit(0);
		}
		System.out.println("Opened and read file: " + filenameFragments);
	}/*  end of getFragmentsFromFile method  */


	/*  start of getIFsFromFile method  */
	private void getIFsFromFile(){
		//  make sure that the number of points has already been determined
		if(!determinedNumberOfPoints){
			System.out.println("Error: the number of points must be determined before reading the IFs as input.\nSystem.exit(0)");
			System.exit(0);
		}
		
		String line = null;
		int lineCounter = 0;
		int minIndex = 0;
		int maxIndex = 0;
		try{
			BufferedReader br = new BufferedReader(new FileReader(inputPath + filenameIFs));			
			while((line = br.readLine()) != null){
				String[] values = line.split("\t");
				if(values.length == 5){
					//  parse the values once that will be used to save time
					int firstIndex = Integer.parseInt(values[0]);
					int secondIndex = Integer.parseInt(values[1]);

					//  initialize the minimum index and the maximum index on the first round
					if(lineCounter == 0){
						minIndex = firstIndex;
						maxIndex = firstIndex;
					}

					//  check for a new minimum index
					if(firstIndex < minIndex){
						minIndex = firstIndex;
					}
					if(secondIndex < minIndex){
						minIndex = secondIndex;
					}

					//  check for a new maximum index
					if(firstIndex > maxIndex){
						maxIndex = firstIndex;
					}
					if(secondIndex > maxIndex){
						maxIndex = secondIndex;
					}
				}
				else{
					//  each row must have in order: index1    index2    IF    stdev    count
					System.out.println("Error: there are not five (5) entries in each row of the IFs file\n"
							+ "System.exit(0)");
					System.exit(0);
				}
				lineCounter++;
			}
			br.close();
		}catch (IOException e) {
			e.printStackTrace();
			System.out.println("Error: cannot find IFs file with filename \"" + filenameIFs + "\"\n" +
			"System.exit(0)");
			System.exit(0);
		}

		if(maxIndex != numberOfPoints){
			System.out.println("Error: the number of points is not the same in the IF values file (" + maxIndex 
					+ ") and the fragments file (" + numberOfPoints + ")\n"
					+ "System.exit(0)");
			System.exit(0);
		}


		IFvalues = new double[numberOfPoints][numberOfPoints];
		for(int i=0; i<numberOfPoints; i++){
			for(int j=0; j<numberOfPoints; j++){
				IFvalues[i][j] = -1;
			}
		}
		IFvariances = new double[numberOfPoints][numberOfPoints];
		for(int i=0; i<numberOfPoints; i++){
			for(int j=0; j<numberOfPoints; j++){
				IFvariances[i][j] = -1;
			}
		}

		line = null;
		try{
			BufferedReader br = new BufferedReader(new FileReader(inputPath + filenameIFs));			
			while((line = br.readLine()) != null){
				String values[] = line.split("\t");
				if(values.length == 5){
					int firstIndex = Integer.parseInt(values[0]) - minIndex;
					int secondIndex = Integer.parseInt(values[1]) - minIndex;

					//  keep the IFvalues matrix symmetric
					IFvalues[firstIndex][secondIndex] = Double.parseDouble(values[2]);
					IFvalues[secondIndex][firstIndex] = IFvalues[firstIndex][secondIndex];
					averageIF += IFvalues[firstIndex][secondIndex];
					numberOfIFValues++;

					//  the IF input file contains the standard deviation for each of the entries
					//  need to square the value to obtain the variance
					double stdev = Double.parseDouble(values[3]);
					if(stdev <= 0){
						IFvariances[firstIndex][secondIndex] = -1;
						IFvariances[secondIndex][firstIndex] = -1;
					}
					else{
						stdev *= stdev;
						//  keep the IFvariances matrix symmetric
						IFvariances[firstIndex][secondIndex] = stdev;
						IFvariances[secondIndex][firstIndex] = stdev;
						defaultVarianceValue += stdev;
						numberOfVarianceValues++;
					}					
				}
				else{
					System.out.println("Error: there are not five (5) entries in each row of the IFs file\n"
							+ "System.exit(0)");
					System.exit(0);
				}
			}
			br.close();
		}catch (IOException e) {
			e.printStackTrace();
			System.out.println("Error: cannot find IFs file with filename \"" + filenameIFs + "\"\n" +
			"System.exit(0)");
			System.exit(0);
		}
		System.out.println("Opened and read file: " + filenameIFs);

		
		//  calculate the average IF value
		if(numberOfIFValues == 0){
			System.out.println("Error: cannot divide by zero. numberOfIFValues is zero (0).\nSystem.exit(0)");
			System.exit(0);
		}
		averageIF = averageIF / numberOfIFValues;
		
		
		//  calculate and set the average variance value for the IF values that do not have a variance value in the input file
		if(numberOfVarianceValues == 0){
			System.out.println("Error: cannot divide by zero. numberOfVarianceValues is zero (0).\nSystem.exit(0)");
			System.exit(0);
		}
		defaultVarianceValue = defaultVarianceValue / numberOfVarianceValues;
		for(int i=0; i<numberOfPoints; i++){
			for(int j=0; j<numberOfPoints; j++){
				if((IFvalues[i][j] != -1) && (IFvariances[i][j] == -1)){
					IFvariances[i][j] = defaultVarianceValue;
				}
			}
		}
	}/*  end of getIFsFromFile method  */


	/*  start of initializeStructure method  */
	private void initializeStructure(){
		
		
		//		CUBE WITH SIDE LENGTH [0,0.01]
		structure = new double[numberOfPoints][3];

		Random random = new Random();
		
		//  initialize the structure in a cube with side length [0,0.01]
		for(int i=0; i<structure.length; i++){
			structure[i][0] = random.nextDouble() / 100.0;
			structure[i][1] = random.nextDouble() / 100.0;
			structure[i][2] = random.nextDouble() / 100.0;
		}
		System.out.println("Initializing the structure in a cube with side length [0,0.01]");
		
		
		
		
		/*
		//		SPRIAL WITH SMALL RANDOM NOISE
		structure = new double[numberOfPoints][3];

		Random random = new Random();
		
		//  use counter1 for the x-coordinates 
		int counter1 = 1;
		//  use counte2 for the y-coordinates
		int counter2 = 2;
		//  use zStepSize for the z-coordinates
		double zStepSize = 5.0 / ((double)numberOfPoints);
		
		//  initialize the structure along a spiral with a bit of random noise
		for(int i=0; i<structure.length; i++){
			//  x-coordinate
			if(counter1 == 1 || counter1 == 2){
				structure[i][0] = 10 + random.nextDouble();
				counter1++;
			}
			else{
				structure[i][0] = 11 + random.nextDouble();
				counter1++;
			}
			//  reset the counter1
			if(counter1 == 5){
				counter1 = 1;
			}
			
			//  y-coordinate
			if(counter2 == 1 || counter2 == 2){
				structure[i][1] = 10 + random.nextDouble();
				counter2++;
			}
			else{
				structure[i][1] = 11 + random.nextDouble();
				counter2++;
			}
			//  reset the counter2
			if(counter2 == 5){
				counter2 = 1;
			}
			
			//  z-coordinate
			structure[i][2] = (i * zStepSize) + random.nextDouble();			
		}
		System.out.println("Initializing the structure in a spiral (with small random noise)");
		*/
		


		
		/*
		//		CUBE WITH SIDE LENGTH = 10 * AVERAGE_DISTANCE
		structure = new double[numberOfPoints][3];

		Random random = new Random();
		
		double averageDistanceX10 = 10.0/averageIF;
		
		//  initialize the structure in a cube with side length = 10 * average distance of input dataset
		for(int i=0; i<structure.length; i++){
			structure[i][0] = random.nextDouble() * averageDistanceX10;
			structure[i][1] = random.nextDouble() * averageDistanceX10;
			structure[i][2] = random.nextDouble() * averageDistanceX10;
			
		}
		System.out.println("Initializing the structure in a cube with side length = " + averageDistanceX10);
		*/
	}/*  end of initializeStructure method  */

	
	/*  start of outputGenomicPositionMap method  */
	private void outputGenomicPositionMap(){
		
		String outputFilename = outputPath + "genomicPositionMap.txt";

		try{						
			BufferedWriter bw = new BufferedWriter(new FileWriter(outputFilename));

			//  write the map
			for(int i=0; i<structure.length; i++){
				bw.write((i+1) + "\t" + pointGenomicLocations[i]);
				bw.newLine();
				bw.flush();				
			}
			bw.flush();
			bw.close();
		}catch (IOException e) {
			System.out.println("Error: could not output file: " + outputFilename);
			e.printStackTrace();
		}//end of try-catch block
		System.out.println("Output file: " + outputFilename);
	}/*  end of outputGenomicPositionMap method  */

	
	/*  start of outputStructure method  */
	private void outputStructure(double structureIndexDouble, double probability, String name){		
		//  define the structure filename
		String outputFilename = new String();
		
		String structureIndexString = formatter.format(structureIndexDouble);
		
		if(name.compareTo("default") == 0){
			outputFilename += "structure_";
			
			//  determine the number of leading zeros
			if(structureIndexDouble < 10){
				outputFilename += "000";
			}
			else if(structureIndexDouble < 100){
				outputFilename += "00";
			}
			else if(structureIndexDouble < 1000){
				outputFilename += "0";
			}
			else{
				//  don't need to do anything
			}
			
			//  append the filetype extension
			outputFilename += structureIndexString + ".pdb";
		}
		else{
			//  use the name passed as a parameter to the method
			outputFilename += name;
		}
		
		
		//  create the header and write the structure to file
		String header = dateFormat.format(new Date()) + "\tlengthOfRun=\t" + formatter.format(runLength) + "\tstepSize=\t" + stepSize +
			"\tindex=\t" + structureIndexString + "\tln(probability)=\t" + probability;
		pdb.writeStructure(outputPath + outputFilename, structure, header);
	}/*  end of outputStructure method  */
	
	
	/*  start of runMCMC method  */
	private void runMCMC(){

		//  initialize the structure
		initializeStructure();
				

		//  output the genomic position map
		outputGenomicPositionMap();

		
		//  set the half step size
		double halfStepSize = stepSize / 2.0;

		
		//  create the random number generator
		Random random = new Random();

		
		//  do not impose any prior biases on the posterior probability of the structures
		double probabilityOfStructure = 1.0;
		double probabilityOfIFData = 1.0;


		//  keep track of the best structure ever seen
		bestStructure = new double[numberOfPoints][3];
		probabilityOfBestStructure = -1;
		indexOfBestStructure = -1;

		try{
			BufferedWriter bw = new BufferedWriter(new FileWriter(outputPath + "stats_" + filenameIFs));

			//  write the header line
			date = new Date();
			String dateString = dateFormat.format(date);
			int mid = filenameIFs.lastIndexOf(".");
			String filenameWithoutExtension = filenameIFs.substring(0,mid);
			bw.write(dateString + "\tlengthOfRun=\t" + formatter.format(runLength) + "\tstepSize=\t" + stepSize + "\tacceptanceFunctionType=\t" + acceptanceFunctionType + "\t" + filenameWithoutExtension);
			bw.flush();
			bw.newLine();
			bw.write("command line args=\t");
			bw.flush();
			for(int i=0; i<argsGlobal.length; i++){
				if(i != argsGlobal.length-1){
					bw.write(argsGlobal[i] + "\t");
				}
				else{
					bw.write(argsGlobal[i]);
				}
				bw.flush();
			}
			bw.newLine();

			
			//  calculate the probability of the initial structure
			//  ln(Pr[S|IF]) = ln(Pr[IF|S]) + ln(Pr[S]) - ln(Pr[IF])
			//	ln(Pr[IF|S]) = SUM(i,j)[- ((IF(i,j) - (1/d(i,j)))^2 / variance)]
			
			double probabilityCurrent = 0;
			
			for(int i=0; i<numberOfPoints-1; i++){
				for(int j=i+1; j<numberOfPoints; j++){
					if(IFvalues[i][j] != -1.0){

						//  calculate the distance in the structure between point i and point j
						double distanceTerm1 = structure[j][0] - structure[i][0];
						distanceTerm1 *= distanceTerm1;
						double distanceTerm2 = structure[j][1] - structure[i][1];
						distanceTerm2 *= distanceTerm2;
						double distanceTerm3 = structure[j][2] - structure[i][2];
						distanceTerm3 *= distanceTerm3;
						
						double distance = distanceTerm1 + distanceTerm2 + distanceTerm3;
						distance = Math.sqrt(distance);
						
						double IFmodel = 1.0 / distance;

						//  calculate the two terms and add the product to the sum
						double term = IFvalues[i][j] - IFmodel;	//  the numerator
						term *= term;
						term /= IFvariances[i][j];		//  divide the numerator by the denominator
						term += Math.log(2.0 * Math.PI * IFvariances[i][j]);

						probabilityCurrent += term;
					}
				}
			}
			probabilityCurrent *= -0.5;
			
			//  output the initial structure
			outputStructure(0, probabilityCurrent, "structure_0000.pdb");
			bw.write(0 + "\t" + probabilityCurrent);
			bw.newLine();
			bw.flush();

			
			//  set the initial structure as the best structure seen
			probabilityOfBestStructure = probabilityCurrent;
			indexOfBestStructure = 0;
			for(int i=0; i<structure.length; i++){
				bestStructure[i][0] = structure[i][0];
				bestStructure[i][1] = structure[i][1];
				bestStructure[i][2] = structure[i][2];
			}
			

			//  use a separate counter so the structures that are output have consecutive numbers in their filenames
			int structureFilenameCounter = 1;
			
			
			//  run the Markov chain Monte Carlo simulation using the Metropolis-Hastings algorithm
			for(long counter=1; counter<=runLength; counter++){

				//  choose a point to move at random
				int pointToMove = random.nextInt(numberOfPoints);
				double currentX = structure[pointToMove][0];
				double currentY = structure[pointToMove][1];
				double currentZ = structure[pointToMove][2];

				
				//  move the point a random distance
				double newX = currentX + (random.nextDouble() * stepSize) - halfStepSize;
				double newY = currentY + (random.nextDouble() * stepSize) - halfStepSize;
				double newZ = currentZ + (random.nextDouble() * stepSize) - halfStepSize;
				

				//  calculate the probability of the new sampled structure
				//  ln(Pr[S|IF]) = ln(Pr[IF|S]) + ln(Pr[S]) - ln(Pr[IF])
				//	ln(Pr[IF|S]) = SUM(i,j)[- ((IF(i,j) - (1/d(i,j)))^2 / variance)]
				double probabilityNew = probabilityCurrent;
								
				for(int j=0; j<numberOfPoints; j++){
					if(j != pointToMove){
						if(IFvalues[pointToMove][j] != -1.0){
							double xj = structure[j][0];
							double yj = structure[j][1];
							double zj = structure[j][2];
							
							//  calculate the distance in the current structure between the point to move (current location) and point j
							double distanceCurrentTerm1 = xj - currentX;
							distanceCurrentTerm1 *= distanceCurrentTerm1;
							double distanceCurrentTerm2 = yj - currentY;
							distanceCurrentTerm2 *= distanceCurrentTerm2;
							double distanceCurrentTerm3 = zj - currentZ;
							distanceCurrentTerm3 *= distanceCurrentTerm3;

							double distanceCurrent = distanceCurrentTerm1 + distanceCurrentTerm2 + distanceCurrentTerm3;
							distanceCurrent = Math.sqrt(distanceCurrent);
							
							double IFmodelCurrent = 1.0 / distanceCurrent;

							//  calculate the distance in the new structure between the point to move (new location) and point j
							double distanceNewTerm1 = xj - newX;
							distanceNewTerm1 *= distanceNewTerm1;
							double distanceNewTerm2 = yj - newY;
							distanceNewTerm2 *= distanceNewTerm2;
							double distanceNewTerm3 = zj - newZ;
							distanceNewTerm3 *= distanceNewTerm3;

							double distanceNew = distanceNewTerm1 + distanceNewTerm2 + distanceNewTerm3;
							distanceNew = Math.sqrt(distanceNew);
							
							double IFmodelNew = 1.0 / distanceNew;

							
							//  make sure not to divide by zero
							if((distanceCurrent <= 0) || (distanceNew <= 0)){
								System.out.println("Error: the distance between the points was not properly determined\n"
										+ "System.exit(0)");
								System.exit(0);
							}

							
							//  update the probability of the new structure (only change the terms that contain the point to move)
							double factor1 = 1.0 / (2 * IFvariances[pointToMove][j]);
							double factor2 = IFmodelCurrent - IFmodelNew;
							double factor3 = IFmodelNew + IFmodelCurrent - (2 * IFvalues[pointToMove][j]);
							double factorUpdate = factor1 * factor2 * factor3;
							
							probabilityNew += factorUpdate;
							
							
							//  slower method
							//  update the probability of the new structure (only change the terms that contain the point to move)							
							//  first subtract the value of the term corresponding to the current location of the point
//							probabilityNew -= (-1.0 / (2.0 * IFvariances[pointToMove][j])) * Math.pow(IFmodelCurrent - IFvalues[pointToMove][j], 2);
							//  second add the value of the term corresponding to the new location of the point
//							probabilityNew += (-1.0 / (2.0 * IFvariances[pointToMove][j])) * Math.pow(IFmodelNew - IFvalues[pointToMove][j], 2);							
						}
					}
				}


				//  accept or reject the new structure using the acceptance function
				double lnProbabilityRatio = probabilityNew - probabilityCurrent;
				double probabilityRatio = Math.pow(Math.E, lnProbabilityRatio);
				double alpha = 1.0;
				if(probabilityRatio < alpha){		//  if the ratio is less than one (1), then the new structure has a lower probability than the current structure
					alpha = probabilityRatio;
					if(acceptanceFunctionGradient){		//  if using the gradient acceptance function, then do NOT want to accept a new structure that is worse than the current structure
						alpha = -1.0;					//  an alpha of negative one (-1) will never be bigger than a randomly chosen positive value and will therefore never be accepted
					}
				}

				if(alpha >= random.nextDouble()){				//  accept the new structure if alpha is bigger than a uniformly chosen random variable [0,1[
					//  update the probability of the structure					
					probabilityCurrent = probabilityNew;
					
					
					//  keep track of the number of changes accepted during the MCMC run
					numberOfChangesAccepted++;


					//  accept the new structure by updating the coordinates of the point to move
					structure[pointToMove][0] = newX;
					structure[pointToMove][1] = newY;
					structure[pointToMove][2] = newZ;

					
					//  check if the new structure is the best structure ever seen
					if(probabilityCurrent > probabilityOfBestStructure){
						probabilityOfBestStructure = probabilityCurrent;
						indexOfBestStructure = counter;

						//  update the coordinates of the best structure
						for(int i=0; i<structure.length; i++){
							bestStructure[i][0] = structure[i][0];
							bestStructure[i][1] = structure[i][1];
							bestStructure[i][2] = structure[i][2];
						}
					}
				}

				//  output the structure and statistics to file
				if((counter % (runLength / numberOfStructuresToOutput)) == 0){				
					outputStructure(structureFilenameCounter, probabilityCurrent, "default");
					bw.write(structureFilenameCounter + "\t" + probabilityCurrent);
					bw.flush();
					bw.newLine();
					structureFilenameCounter++;
				
					//  output and overwrite the best structure seen so far
					double[][] tempStructure = structure;
					structure = bestStructure;				//  note: the outputStructure method will output the structure pointed to by the structure global variable
					outputStructure(indexOfBestStructure, probabilityOfBestStructure, "best_structure.pdb");
					structure = tempStructure;					
				}
			}//end of Markov chain Monte Carlo loop
			bw.close();
		}catch (IOException e) {
			System.out.println("Error: could not output file: " + outputPath + "stats_" + filenameIFs);
			e.printStackTrace();
		}//end of try-catch block
		System.out.println("Output file: " + outputPath + "stats_" + filenameIFs);
		
		
		//  output the best structure seen during the run
		structure = bestStructure;				//  note: the outputStructure method will output the structure pointed to by the structure global variable
		outputStructure(indexOfBestStructure, probabilityOfBestStructure, "best_structure.pdb");
	}/*  end of runMCMC method  */


	/*  start of run method  */
	public void run(String[] args) {
		argsGlobal = args;
		if(args.length == 2){
			filenameIFs = args[0];
			filenameFragments = args[1];
		}
		else if(args.length == 3){
			filenameIFs = args[0];
			filenameFragments = args[1];
			runLength = Long.parseLong(args[2]);
		}
		else if(args.length == 4){
			filenameIFs = args[0];
			filenameFragments = args[1];
			runLength = Long.parseLong(args[2]);
			numberOfStructuresToOutput = Double.parseDouble(args[3]);
		}
		else if(args.length == 5){
			filenameIFs = args[0];
			filenameFragments = args[1];
			runLength = Long.parseLong(args[2]);
			numberOfStructuresToOutput = Double.parseDouble(args[3]);
			stepSize = Double.parseDouble(args[4]);
		}
		else if(args.length == 6){
			filenameIFs = args[0];
			filenameFragments = args[1];
			runLength = Long.parseLong(args[2]);
			numberOfStructuresToOutput = Double.parseDouble(args[3]);
			stepSize = Double.parseDouble(args[4]);
			acceptanceFunctionType = args[5];
			if((acceptanceFunctionType.compareTo("mcmc") != 0) && (acceptanceFunctionType.compareTo("gradient") != 0)){
				System.out.println("Error: wrong type of acceptance function specified. Use \"mcmc\" or \"gradient\"\nSystem.exit(0)");
				System.exit(0);
			}
		}
		else{
//			System.out.println("Error:\tplease use one of the following patterns (without the quotes)\n" +
//					"\t\"java MCMC5CMain IFsFilename.txt fragmentsFilename.txt\"\n" +
//					"\t\"java MCMC5CMain IFsFilename.txt fragmentsFilename.txt runLength\"\n" +
//					"\t\"java MCMC5CMain IFsFilename.txt fragmentsFilename.txt runLength numberOfStructuresToOutput\"\n" +
//					"\t\"java MCMC5CMain IFsFilename.txt fragmentsFilename.txt runLength numberOfStructuresToOutput stepSize\"\n" +
//					"\t\"java MCMC5CMain IFsFilename.txt fragmentsFilename.txt runLength numberOfStructuresToOutput stepSize acceptanceFunctionType\"\n" +
//				"Notes:\tthe default run length is " + formatter.format(runLength) + "\n" +
//					"\tthe default number of structures to output is " + (int)numberOfStructuresToOutput + "\n" +
//					"\tthe default variance of the noise model is " + defaultVarianceValue + "\n" +
//					"\tthe default step size is " + stepSize + "\n" +
//					"\tthe default acceptance function type is " + acceptanceFunctionType + "\n" +
//				"System.exit(0)");
//			System.exit(0);
		}
		
		
		
		
		
		
		
		
		
		/** ------------------------------------------------------------------------------------------------------------------
		 * 
		 * 
		 * 								HARDCODING OF PARAMETERS
		 * 
		 * 
		 * ------------------------------------------------------------------------------------------------------------------  */
		filenameIFs = "DMSO_experimental.txt";
//		filenameIFs = "PMA_experimental.txt";
		filenameFragments = "fragments.txt";

//		filenameIFs = "HiC_chr14_chr14.txt";
//		filenameFragments = "HiC_fragments.txt";

		runLength = 10000000;

		numberOfStructuresToOutput = 10;

		stepSize = 0.05;

		acceptanceFunctionType = "mcmc";
//		acceptanceFunctionType = "gradient";
		/** ------------------------------------------------------------------------------------------------------------------
		 * 
		 * 
		 * 								END OF HARDCODING OF PARAMETERS
		 * 
		 * 
		 * ------------------------------------------------------------------------------------------------------------------  */
				
		
		
		
		


				
		//  create the output folder
		new File(outputFolder).mkdir();
		
		
		
		//  output the parameters to file
		try{
			String parameterFilename = "parameters.txt";
			BufferedWriter bw = new BufferedWriter(new FileWriter(outputPath + parameterFilename));
			
			bw.write("filenameIFs=\t" + filenameIFs);
			bw.newLine();
			bw.flush();
			bw.write("filenameFragments=\t" + filenameFragments);
			bw.newLine();
			bw.flush();
			bw.write("runLength=\t" + runLength);
			bw.newLine();
			bw.flush();
			bw.write("numberOfStructuresToOutput=\t" + (int)numberOfStructuresToOutput);
			bw.newLine();
			bw.flush();
			bw.write("stepSize=\t" + stepSize);
			bw.newLine();
			bw.flush();
			bw.write("acceptanceFunctionType=\t" + acceptanceFunctionType);
			bw.newLine();
			bw.flush();
			bw.write("command line parameters:");
			for(String arg: args){
				bw.write("\t" + arg);
			}
			bw.newLine();
			bw.flush();
			
			bw.close();
			System.out.println("Output file: " + outputPath + parameterFilename);
		}catch (IOException e) {
			e.printStackTrace();
		}//end of try-catch block
		

		
		//  get the number of points (can pass the fragments file as parameter to count the number of points from this file)
		//  the number of points MUST be determined before the fragments or the IFs are read as input
		getNumberOfPoints(filenameFragments);

		
		//  get the fragment genomic locations
		getFragmentsFromFile();


		//  get the IF values (must be done after getting the fragment genomic locations)
		getIFsFromFile();


		//  output the parameters of the run
		System.out.println("Run parameters: ");
		System.out.println("\tLength of run = " + runLength);
		System.out.println("\tNumber of structures to output = " + (int)numberOfStructuresToOutput);
		System.out.println("\tThe default variance of the noise model = " + defaultVarianceValue);
		System.out.println("\tStep size = " + stepSize);
		System.out.println("\tAcceptance function type = " + acceptanceFunctionType);
		
		if(acceptanceFunctionType.compareTo("mcmc") == 0){
			acceptanceFunctionMCMC = true;
			acceptanceFunctionGradient = false;
		}
		else if(acceptanceFunctionType.compareTo("gradient") == 0){
			acceptanceFunctionMCMC = false;
			acceptanceFunctionGradient = true;
		}
		

		//  calculate how long running the MCMC takes (time point before run)
		long timeBefore = System.currentTimeMillis();


		//  run the MCMC simulation
		runMCMC();


		//  output the time taken by the MCMC run
		long timeTaken = System.currentTimeMillis() - timeBefore;
		timeTaken = timeTaken / (long)(1000.0);			//	convert from milliseconds to seconds		
		if(timeTaken > 3600){		//	(1hr = 3600 sec)
			int numberOfHours = (int)timeTaken / (int)3600;
			long leftOver = timeTaken % 3600;
			int numberOfMinutes = (int)leftOver / (int)60;
			leftOver = leftOver % 60;
			int numberOfSeconds = (int)leftOver;
			System.out.println("Amount of time taken = " + numberOfHours + "hr" + numberOfMinutes + "min" + numberOfSeconds + "sec");
		}
		else if(timeTaken > 60){		//	(1min = 60 sec)
			int numberOfMinutes = (int)timeTaken / (int)60;
			long leftOver = timeTaken % 60;
			int numberOfSeconds = (int)leftOver;
			System.out.println("Amount of time taken = " + numberOfMinutes + "min" + numberOfSeconds + "sec");
		}
		else{
			int numberOfSeconds = (int)timeTaken;
			System.out.println("Amount of time taken = " + numberOfSeconds + "sec");
		}

		
		//  output the number of proposed changes accepted during the run
		System.out.println("Number of proposed changes accepted: " + numberOfChangesAccepted +
				" [" + (int)(numberOfChangesAccepted * 100 / runLength) + "%]");
		
		System.out.println("Done!");
	}/*  end of run method  */
}/*  end of MCMC5C class  */